# -*- coding: utf-8 -*-
"""
Created on Tue Feb  5 13:35:30 2013

@author: daniel
"""
import numpy as np
from scipy.optimize import fsolve
from scipy.optimize import minimize
class PBDO:
    def __init__(self,obj,ineqcons,params,args=(),alpha=0.0,dx=10**-7,nCPU=1,maxIter=100,tol=10**-5,iprint=1,optimizer='SLSQP'):
        self.optimizer=optimizer        
        """
        PBDO Optimization class
        
        Parameters
        ----------
        
        obj : function
            objective function
        ineqcons : list of functions
            inequality constraints list
        params : parameters
            parameters class
        args : list
            additional arguments to be passed to objective function
        alpha : float
            alpha value
        dx : float
            dx value for gradient determination
        nCpu : integer
            number of CPU for parallel computation
        maxIter : integer
            maximum number of iterations
        tol : float
            tolerance
        iprint : integer
            display the result
        
        """
        self.nCPU=1 #Temporarily restricted to 1 CPU core        
        self.dx=dx  
        self.tol=tol
        self.obj=obj        
        self.args=args
        self.cons=ineqcons        
        
        self.alpha=alpha
        self.params=params
        
        self.d0=self.params.deterministicVars.d0
        self.dl=self.params.deterministicVars.lb
        self.du=self.params.deterministicVars.ub
        
        self.x0=self.params.fuzzyVars.x0
        self.xl=self.params.fuzzyVars.lb
        self.xu=self.params.fuzzyVars.ub
        
        self.P0=self.params.fuzzyParams.x0       
        self.S =list()
        self.Pmpp=list()
        self.xHistory=list()
        self.fHistory=list()
        self.dHistory=list()
        self.dMin=np.zeros(len(self.d0))
        self.xMin=np.zeros(len(self.x0))
        self.fMin=0.0
        self.xHistory_full=list()
        self.fHistory_full=np.array([])
        for i in range(len(self.cons)):
            self.S.append(np.zeros([len(self.x0)]))
            self.Pmpp.append(self.P0)
    def solve(self):
        err=np.inf
        iMax=5
        i=1
        while err>self.tol and i<=iMax:
            f1=self.obj(self.d0,self.x0,self.P0)
            self.deterministicPhase()
            self.possibilityAssessmentPhase()
            f2=self.obj(self.d0,self.x0,self.P0)
            err=np.abs(f2-f1)
            i+=1
            self.xHistory.append(self.x0)
            self.fHistory.append(f2)
            self.dHistory.append(self.d0)
    def set_d0(self,d0):
        assert len(d0)==len(self.params.deterministicVars.d0)
        assert len(d0)==len(self.params.deterministicVars.lb)
        assert len(d0)==len(self.params.deterministicVars.ub)
        self.params.deterministicVars.d0=d0
    def set_x0(self,x0):
        assert len(self.fuzzyVars)==len(x0)
        for i in range(len(x0)):self.params.fuzzyVars[i].set_x0(x0)
    def set_alpha(self,alpha):
        self.alpha=alpha
    def deterministicPhase(self):
        """
        Runs a deterministic optimizaiton using shifted constraint and 
        parameter values based on the possibility analysis of each constraint
        """
        nd=len(self.d0)  # Number of deterministic design variables
        nx=len(self.x0)  # Number of fuzzy design variables
        
        dx0 =np.concatenate((self.d0,self.x0)) # joining fuzzy and deterministic
        dxLB=np.concatenate((self.dl,self.xl)) # variables and boundaries
        dxUB=np.concatenate((self.du,self.xu))
        
        mf=self.modifiedFunction(self.obj,nd,nx,np.zeros(len(self.x0)),self.P0)   
        # Setting the modified objective funciton and constraint functions so
        # they are compatable with scipy optimizers
        obj=mf.modFcn
        cons=list()
        
        for i,con in enumerate(self.cons):
            mf=self.modifiedFunction(con,nd,nx,self.S[i],self.Pmpp[i])
            conFcn=mf.modFcn
            cons.append({'type':'ineq','fun':conFcn})
        bounds=self._getBounds(dxLB,dxUB)
        # Calling a multicore wrapper for scipy.slsqp
        #fmin=fmin_multicore.fmin(obj,dx0,bounds,dx=self.dx,ineqcons=cons,args=self.args,nCPU=self.nCPU,iprint=0)        
        fmin         =minimize(obj,dx0,bounds=bounds,constraints=cons,args=self.args,method=self.optimizer)
        #fmin.solve()
        dxMin=fmin.x
        d=dxMin[0:nd]
        x=dxMin[nd:]
        self.d0=d
        self.x0=x
        self.xMin=x
        self.dMin=d
        self.fMin=fmin.fun
#        self.xHistory_full.append(fmin.xHistory)
#        self.fHistory_full=np.concatenate((self.fHistory_full,fmin.fHistory))
    def _getBounds(self,lb,ub):
        assert len(lb)==len(ub)
        bounds=np.zeros([len(lb),2])
        for i in range(len(lb)):
            bounds[i,0]=lb[i]
            bounds[i,1]=ub[i]
        return bounds
    class modifiedFunction:
        def __init__(self,fcn,nd,nx,S,P):
            self.fcn=fcn    
            self.nd=nd
            self.nx=nx
            self.P=P
            self.S=S
        def modFcn(self,dx,*args):
            d=dx[0:self.nd]
            x=dx[self.nd:]
            return self.fcn(d,x+self.S,self.P,*args)
    def possibilityAssessmentPhase(self):
        nx=len(self.x0)
        nP=len(self.P0)
        xlb,x0,xub=self.params.fuzzyVars.getVectors(self.alpha,self.x0)
        Plb,P0,Pub=self.params.fuzzyParams.getVectors(self.alpha)
        
        xP0 =np.concatenate((x0,P0))
        xPLB=np.concatenate((xlb,Plb))
        xPUB=np.concatenate((xub,Pub))
        bounds=self._getBounds(xPLB,xPUB)
        for i,con in enumerate(self.cons):
            PMAfcn       =self.PMA_modifiedFunctions(con,self.d0,nx,nP)
            PMAobj       =PMAfcn.PMAobj
            #fmin         =fmin_multicore.fmin(PMAobj,xP0,bounds,dx=self.dx,args=self.args,nCPU=self.nCPU,iprint=0)
            fmin         =minimize(PMAobj,xP0,bounds=bounds,args=self.args,method=self.optimizer)
#            fmin.solve()
            xPmin        =fmin.x
            xMPP         =xPmin[0:nx]
            Pmpp         =xPmin[nx:]
            self.Pmpp[i] =Pmpp
            self.S[i]    =xMPP-self.x0
    class PMA_modifiedFunctions:
        def __init__(self,con,d,nx,nP):
            self.con     =con
            self.nx      =nx
            self.nP      =nP
            self.d       =d
        def PMAobj(self,xP,*args):
            x           =xP[0:self.nx]
            P           =xP[0:self.nP]
            f           =self.con(self.d,x,P)
            return f
    def display(self):
        for i,f in enumerate(self.fHistory):
            print '{:<4d}f = {:<10.5f}  '.format(i,f),
            print 'x=[',
            for j,x in enumerate(self.xHistory[i]):
                print '{:<6.4f}'.format(x),
            print '] d=[',
            for k,d in enumerate(self.dHistory[i]):
                print '{:<6.4f}'.format(d),
            print ']'
        self.fHistory=list()
        self.xHistory=list()
        self.dHistory=list()
         
class parameters:
    def __init__(self):
        self.fuzzyParams=fuzzyVector()
        self.fuzzyVars=fuzzyVector()
        self.deterministicVars=deterministicVars()
    def addDeterminisitcVar(self,d,lb,ub):
        self.deterministicVars.addVariable(d,lb,ub)
    def addFuzzyVar(self,fuzzyVar,lb,ub):
        self.fuzzyVars.add(fuzzyVar,lb,ub)
    def addFuzzyParam(self,fuzzyParam):
        self.fuzzyParams.add(fuzzyParam)
    
class deterministicVars:
    def __init__(self):
        self.d0=np.array([])
        self.lb=np.array([])
        self.ub=np.array([])
    def addVariable(self,d,lb,ub):
        assert len(d)==len(lb)
        assert len(lb)==len(ub)
        if type(d) ==type(float()):d =np.array([d])
        if type(lb)==type(float()):lb=np.array([lb])
        if type(ub)==type(float()):ub=np.array([ub])
        self.d0=np.concatenate((self.d0,d))
        self.lb=np.concatenate((self.lb,lb))
        self.ub=np.concatenate((self.ub,ub))
    def __getitem__(self,k):
        return (self.d0[k],self.lb[k],self.ub[k])
              
class fuzzyVector:
    def __init__(self):
        self._list=list()
        self.lb=np.array([])
        self.ub=np.array([])
        self.x0=np.array([])
    def add(self,fuzzyNumber,lb=None,ub=None):
        self._list.append(fuzzyNumber)
        self.lb=np.concatenate((self.lb,np.array([lb])))
        self.ub=np.concatenate((self.ub,np.array([ub])))
        self.x0=np.concatenate((self.x0,np.array([fuzzyNumber.x0])))
    def getVectors(self,alpha,x=None):
        lb=np.zeros(len(self._list))
        ub=np.zeros(len(self._list))
        X=np.zeros(len(self._list))
        for i,fn in enumerate(self._list):
            if x==None:
                I=fn.getInterval(alpha)
            else:
                I=fn.getInterval(alpha,x[i])
            lb[i]=I[0]
            X[i] =I[1]
            ub[i]=I[2]
        return lb,X,ub
    def display(self):
        lb,x,ub=self.getVectors(0.0)
        print lb
        print x
        print ub
    def __getitem__(self,k):
        return self._list[k]


class fuzzyNumber:
    """
    Class containing functions for defining and solving fuzzy numbers
    The input function must comply with the follwing definitions:
    1. Must have a maximal value of 1 at some x, called x0
    2. Must have an x-intercept at x<x0 called xMin
    3. Must have an x-intercept at x>x0 called xMax
    4. Must be continuous between xMin and xMax
    """
    def __init__(self,fcn):
        self.function=fcn
        self._delta=0.0
        self._dx=0.001
        self._x0=0.0
        self.xMin=0.0
        self.xMax=0.0
        self._delta=0.0
        self._find_x0()
        self._find_xMin()
        self._find_xMax()
        self.x0=self._x0
    def set_x0(self,x0):
        self.x0=x0
        self._delta=x0-self._x0
    def _find_x0(self):
        self._x0=fsolve(self._x0Eqn,0.0)[0]
    def _find_xMin(self):
        self.xMin=fsolve(self.function,self._x0-self._dx)[0]
    def _find_xMax(self):
        self.xMax=fsolve(self.function,self._x0+self._dx)[0]
    def _x0Eqn(self,x):
        return self.function(x)-1.0
    def getInterval(self,alpha,x0=None):
        if x0==None:x0=self.x0
        else:self.set_x0(x0)
        dxR=(self.xMax+self._delta-self.x0)/2.
        dxL=(self.x0-self.xMin-self._delta)/2.
        if alpha==0:
            interval=np.array([self.xMin+self._delta,x0,self.xMax+self._delta])
        elif alpha <1.0 and alpha>0.0:
            xMin0=fsolve(self._xEqn,self.xMin-dxL,args=alpha)[0]
            xMax0=fsolve(self._xEqn,self.xMax+dxR,args=alpha)[0]
            interval=np.array([xMin0+self._delta,x0,xMax0+self._delta])
        elif alpha==1.0:
            interval=np.array([x0,x0,x0])
        else:
            interval=None
            print 'Error, 0<=alpha<=1'
        return interval
    def _xEqn(self,x,alpha):
        return self.function(x)-alpha
    def getV(self,x):
        I=self.getInterval(0.0)
        V=(x-I[0])/(I[1]-I[0])*2.0-1.0
        return V
    def getAlpha(self,x):
       x0=x-self._delta
       alpha=self.function(x0)
       return alpha


class triangle:
    """
    Creates a triangular fuzzy function as follows
              _______________alpha=1
             /|\
            / | \ 
           /  |  \
          /   |   \
     lb__/    x    \__ub_ alhpa=0         
    """
    def __init__(self):
        self.lb=0.0
        self.ub=0.0
        self.x =0.0
    def getFuzzyNumber(self,lb,x0,ub):
        self.lb=lb
        self.ub=ub
        self.x=x0
        return fuzzyNumber(self.fcn)
    def fcn(self,x):
        mLHS=+1.0/(self.x-self.lb)
        bLHS=-self.lb*mLHS
        mRHS=-1.0/(self.ub-self.x)
        bRHS=-self.ub*mRHS
        if x<self.x:
            y=mLHS*x+bLHS
        if x>self.x:
            y=mRHS*x+bRHS
        if x==self.x:
            y=1.0
        return y

class parabola:
    def __init__(self):
        self.lb = 0.0
        self.ub = 0.0
        self.x = 0.0
    def getFuzzyNumber(self,lb,x0,ub):
        self.lb = lb
        self.x= x0
        self.ub = ub
        return fuzzyNumber(self.fcn)
    
    def fcn(self,x):
        pass
#------------------------------------------------------------------------------
def fuzzyNumberTest():
    def fcn(x):
        return 1.0-(x-2.0)**2
    fn=fuzzyNumber(fcn)
    print fn.x0
    print fn.xMin
    print fn.xMax
    print fn.getInterval(0)
    print fn.getInterval(0.4,5.0)
    fn.set_x0(2.0)
    print fn.getInterval(0.2)
    tri=triangle()
    fn=tri.getFuzzyNumber(-1.0,0.0,2.0)

    print fn.getInterval(0.0)
    print fn.getInterval(0.2)
    print fn.getInterval(0.9)
    print fn.getInterval(0.9,2)
    print fn.getInterval(0.9)
    fn.set_x0(5.)
    print fn.getInterval(0.9)
    print fn.getInterval(0)
    print fn.getV(4.1)  
    print fn.getAlpha(4.1)

    tri = triangle()
    fn2 = tri.getFuzzyNumber(0,5,10)
    print fn2.getInterval(0.97)
    
def PBDOtest():
    tri=triangle()    
    params=parameters()
    params.addFuzzyVar(tri.getFuzzyNumber(3.6,5.0,6.4),0,10)
    tri=triangle()    
    params.addFuzzyVar(tri.getFuzzyNumber(4.3,5.0,5.7),0,10)
    opt=PBDO(test1_ob,(test1_C1,test1_C2,test1_C3),params,dx=10**-10,alpha=0.02,nCPU=1)  
    #opt.possibilityAssessmentPhase()
    opt.solve()
    print 'Example 1'
    opt.display()
    
    params=parameters()
    tri=triangle()    
    params.addFuzzyParam(tri.getFuzzyNumber(900.,1000.,1100.))
    tri=triangle()    
    params.addFuzzyParam(tri.getFuzzyNumber(400.,500.,600.))
    tri=triangle()
    params.addFuzzyParam(tri.getFuzzyNumber(38000.,40000.,42000.))
    tri=triangle()
    params.addFuzzyParam(tri.getFuzzyNumber(2.755,2.900,3.045))
    tri=triangle()   
    params.addDeterminisitcVar([1.,3.],[0.1,0.1],[10.,10.]) 
        
    params.addDeterminisitcVar([3.,3.],[0.1,0.1],[10.,10.])

    cons=list()
    cons.append(test2_C1)
    cons.append(test2_C2)


    opt=PBDO(test2_ob,cons,params,dx=10**-7,alpha=0.05,nCPU=1) 
    #opt.possibilityAssessmentPhase()
    opt.solve()
    print 'Example 2'
    opt.display()
#--------------Test Problem 1---------------    
def test1_ob(d,x,P):
    f=x[0]+x[1]
    return f
def test1_C1(d,x,P):
    C1=x[0]**2*x[1]/20.0-1.
    return C1
def test1_C2(d,x,P):
    C2=(x[0]+x[1]-5.0)**2/30.+(x[0]-x[1]-12.0)**2/120.0-1.
    return C2
def test1_C3(d,x,P):
    C3=80.0/(x[0]**2+8.0*x[1]+5.0)-1
    return C3
#--------------Test Problem 2---------------
def test2_ob(d,x,P):
    f=d[0]*d[1]
    return f
def test2_C1(d,x,P):
    C1=P[2]-(600./d[0]/d[1]**2*P[0]+600./d[0]**2/d[1]*P[1])
    return C1
def test2_C2(d,x,P):
    D0=2.5
    L=100.0
    C2=D0-4.0*L**3/(P[3]*10.**7)/d[0]/d[1]*((P[0]/d[1]**2)**2+(P[1]/d[0]**2)**2)**.5
    return C2


if __name__=='__main__':
    fuzzyNumberTest()
    #PBDOtest()
        